#### SETUP #### 

## Packages
need <- c("foreign", "stargazer", "arm", "openxlsx", "clubSandwich", "nnet", "margins", "dplyr", "mediation", "systemfit", "ggplot2", "jtools")
have <- need %in% rownames(installed.packages()) 
if(any(!have)) install.packages(need[!have]) # If this gives you error, this might be your sign that your packages are actually not up-to-date. 
invisible(lapply(need, library, character.only=T))
rm(list = ls()) 

## Opening the datasets
setwd("~/Google Drive File Stream/My Drive/LebFLS_DiversityQuality/Accepted Version") # change the working directory based on where you stored the data
data.tot.iop.leb <- read.xlsx("AllCentersLebanesePatients.xlsx")
data.tot.iops.leb <- read.xlsx("SectarianCentersLebanesePatients.xlsx")
data.tot.iops.syr <- read.xlsx("SectarianCentersSyrianPatients.xlsx")

#### DESCRIPTIVES IN DATASET INTRODUCTION ####

# Descriptives in the dataset introduction part [Table 1]
table(data.tot.iops.leb$centersect, useNA = "always") 




#### SECTION: SELECTION INTO INGROUP AND OUTGROUP CENTERS #### 

## Creating a Sankey diagram [Figure 1]
table(data.tot.iops.leb$sect,data.tot.iops.leb$centersect, useNA = 'always')
# Sankey diagramm was then created using the online tool at http://sankeymatic.com/ 



## For Table 2
tapply(data.tot.iops.leb$inout, data.tot.iops.leb$centersect, sum, na.rm=T)
tapply(data.tot.iops.leb$outin, data.tot.iops.leb$centersect, sum, na.rm=T)
tapply(data.tot.iops.leb$inout, data.tot.iops.leb$centersect, mean, na.rm=T)
mean(data.tot.iops.leb$inout, na.rm=T)
tapply(data.tot.iops.leb$inout, data.tot.iops.leb$centersect, sd, na.rm=T)
sd(data.tot.iops.leb$inout, na.rm=T)
tapply(data.tot.iops.leb$inout, data.tot.iops.leb$centersect, var, na.rm=T)
var(data.tot.iops.leb$inout, na.rm=T)

# for within-center variance
anova(lm(data=data.tot.iops.leb[(data.tot.iops.leb$centersect=="armencenter"),], inout ~ centercode))["Residuals", "Mean Sq"]
anova(lm(data=data.tot.iops.leb[(data.tot.iops.leb$centersect=="christcenter"),], inout ~ centercode))["Residuals", "Mean Sq"]
anova(lm(data=data.tot.iops.leb[(data.tot.iops.leb$centersect=="druzecenter"),], inout ~ centercode))["Residuals", "Mean Sq"]
anova(lm(data=data.tot.iops.leb[(data.tot.iops.leb$centersect=="shiacenter"),], inout ~ centercode))["Residuals", "Mean Sq"]
anova(lm(data=data.tot.iops.leb[(data.tot.iops.leb$centersect=="sunnicenter"),], inout ~ centercode))["Residuals", "Mean Sq"]
anova(lm(data=data.tot.iops.leb, inout ~ centercode))["Residuals", "Mean Sq"]

# for between-center variance
anova(lm(data=data.tot.iops.leb[(data.tot.iops.leb$centersect=="armencenter"),], inout ~ centercode))["centercode", "Mean Sq"]
anova(lm(data=data.tot.iops.leb[(data.tot.iops.leb$centersect=="christcenter"),], inout ~ centercode))["centercode", "Mean Sq"]
anova(lm(data=data.tot.iops.leb[(data.tot.iops.leb$centersect=="druzecenter"),], inout ~ centercode))["centercode", "Mean Sq"]
anova(lm(data=data.tot.iops.leb[(data.tot.iops.leb$centersect=="shiacenter"),], inout ~ centercode))["centercode", "Mean Sq"]
anova(lm(data=data.tot.iops.leb[(data.tot.iops.leb$centersect=="sunnicenter"),], inout ~ centercode))["centercode", "Mean Sq"]
anova(lm(data=data.tot.iops.leb, inout ~ centercode))["centercode", "Mean Sq"]



## For Figure 2
# Panel A
gd <- data.tot.iops.leb %>% 
  group_by(affil) %>% 
  summarise(
    prox = sum(prox, na.rm=TRUE),
    afrd = sum(afrd, na.rm=TRUE),
    cent_rept = sum(cent_rept, na.rm=TRUE),
    doc_rept = sum(doc_rept, na.rm=TRUE),
    doc_know = sum(doc_know, na.rm=TRUE),
    doc_ppl = sum(doc_ppl, na.rm=TRUE),
    specialty = sum(specialty, na.rm=TRUE), 
  )
gd <- gd[1:2,]
# First row of "gd" to be divided by 226 (number of ingroup patients) 
# and second row of "gd" to be divided by 38 (number of outgroup patients). 
table(data.tot.iops.leb$outin)

# Panel B
gd2 <- data.tot.iops.leb[(data.tot.iops.leb$access_io==1),] %>% 
  group_by(affil) %>% 
  summarise(
    prox = sum(prox, na.rm=TRUE),
    afrd = sum(afrd, na.rm=TRUE),
    cent_rept = sum(cent_rept, na.rm=TRUE),
    doc_rept = sum(doc_rept, na.rm=TRUE),
    doc_know = sum(doc_know, na.rm=TRUE),
    doc_ppl = sum(doc_ppl, na.rm=TRUE),
    specialty = sum(specialty, na.rm=TRUE), 
  )
gd2 <- gd2[1:2,]
# Second row of "gd2" to be divided by 151 (number of ingroup patients with access to both ingroup and outgroup centers) 
# and second row of "gd" to be divided by 38 (number of outgroup patients with access to both ingroup and outgroup centers). 
table(data.tot.iops.leb[(data.tot.iops.leb$access_io==1),]$outin)



## Online Appendix Section 2: Descriptive statistics of the center selection variables 
stargazer(data.tot.iops.leb[c("prox", "afrd", "cent_rept", "doc_rept", "doc_know", "doc_ppl", "specialty")],
          type="text",
          covariate.labels = c("Proximity", "Affordability", "Center reputation", "Doctor reputation", 
                               "Knowing medical staff", "Knowing brokers", "Special care"))
stargazer(data.tot.iops.leb[(data.tot.iops.leb$access_io==1),][c("prox", "afrd", "cent_rept", "doc_rept", "doc_know", "doc_ppl", "specialty")],
          type="text",
          covariate.labels = c("Proximity", "Affordability", "Center reputation", "Doctor reputation", 
                               "Knowing medical staff", "Knowing brokers", "Special care"))
stargazer(data.tot.iop.leb[c("prox", "afrd", "cent_rept", "doc_rept", "doc_know", "doc_ppl", "specialty")],
          type="text",
          covariate.labels = c("Proximity", "Affordability", "Center reputation", "Doctor reputation", 
                               "Knowing medical staff", "Knowing brokers", "Special care"))


## Online Appendix Section A3: Multinomial and logistic analyses
# With demographic control variables
test <- multinom(reltype ~ prox + afrd + cent_rept + doc_rept + specialty + doc_know + doc_ppl + 
                   age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious, 
                 data = data.tot.iop.leb)
test.2 <- glm(outin ~ prox + afrd + cent_rept + doc_rept + specialty + doc_know + doc_ppl +
                age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious,
              family=binomial(link="logit"),
              data = data.tot.iops.leb)
test.3 <- glm(outin ~ prox + afrd + cent_rept + doc_rept + specialty + doc_know + doc_ppl +
                age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious,
              family=binomial(link="logit"),
              data = data.tot.iops.leb[(data.tot.iops.leb$access_io==1),])

# Block-bootstrapped standard errors (for multinomial logit)
resamp.block <- function(data){
  lookup <- split(1:nrow(data), data$providercode_unique)
  providers <- names(lookup)
  providers.drawn <- sample(providers, length(unique(data$providercode_unique)), replace = TRUE)
  star <- unlist(lookup[providers.drawn])
  samp.blo.data <- data[star,]
}

coefs.non <- matrix(NA, 15, 1000)
coefs.sect <- matrix(NA, 15, 1000)

for(i in 1:1000){
  resamp.blo.data <- resamp.block(data.tot.iop.leb)
  pref.fit1 <- multinom(reltype ~ prox + afrd + cent_rept + doc_rept + specialty + doc_know + doc_ppl +
                          age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious, 
                        data = resamp.blo.data)
  coefs.non[,i] <- summary(pref.fit1)$coefficients[1,]
  coefs.sect[,i] <- summary(pref.fit1)$coefficients[2,]
  if(i %% 100==0) {
    cat(paste0("iteration: ", i, "\n"))
  }
}

se.pref.fit1a <- list(c(sd(coefs.non[1,]), sd(coefs.non[2,]), 
                        sd(coefs.non[3,]), sd(coefs.non[4,]), 
                        sd(coefs.non[5,]), sd(coefs.non[6,]), 
                        sd(coefs.non[7,]), sd(coefs.non[8,]),
                        sd(coefs.non[9,]), sd(coefs.non[10,]), 
                        sd(coefs.non[11,]), sd(coefs.non[12,]),
                        sd(coefs.non[13,]), sd(coefs.non[14,]), 
                        sd(coefs.non[15,])))
se.pref.fit1b <- list(c(sd(coefs.sect[1,]), sd(coefs.sect[2,]), 
                        sd(coefs.sect[3,]), sd(coefs.sect[4,]), 
                        sd(coefs.sect[5,]), sd(coefs.sect[6,]), 
                        sd(coefs.sect[7,]), sd(coefs.sect[8,]),
                        sd(coefs.sect[9,]), sd(coefs.sect[10,]), 
                        sd(coefs.sect[11,]), sd(coefs.sect[12,]), 
                        sd(coefs.sect[13,]), sd(coefs.sect[14,]), 
                        sd(coefs.sect[15,])))

round(unlist(se.pref.fit1a), 3)
round(unlist(se.pref.fit1b), 3)

# Clustering-robust standard errors for logistic analyses
se.test.2 <- list(coef_test(test.2, vcov = "CR2", 
                            cluster = data.tot.iops.leb$providercode_unique, test = "Satterthwaite")[,2])
se.test.3 <- list(coef_test(test.3, vcov = "CR2", 
                            cluster = data.tot.iops.leb[(data.tot.iops.leb$access_io==1),]$providercode_unique, test = "Satterthwaite")[,2])

stargazer(test, test.2, test.3, 
          type="text",
          covariate.labels = c("Proximity", "Affordability", "Center reputation", "Doc reputation", 
                               "Special care", "Knowing medical staff", "Knowing brokers",
                               "Patient age", "Respondent age", "Patient gender: Female", "Respondent gender: Female", 
                               "Socioecon. status", "Patient general health", "Religiosity"),
          column.labels = c("All Lebanese patients" , "All Lebanese patients", "Lebanese patients in sectarian centers" , "Lebanese patients in sectarian centers with access to both ingroup and outgroup centers"),
          dep.var.labels = c("Non-sectarian center vs. ingroup center", "Sectarian outgroup center vs. ingroup center", 
                             "Outgroup center vs. ingroup center", "Outgroup center vs. ingroup center"),
          notes.align = "l",
          notes.append = FALSE,
          notes = c("* denotes p-value < 0.1, ** denotes p-value < 0.05 and *** denotes p-value < 0.01.", 
                    "Normal standard errors are presented in parentheses.",
                    "Cluster-robust standard errors (at the provider level) are presented in brackets."))
# Manual changes in the table: Enter the cluster-robust standard errors manually. 



## Online Appendix Section 4: SUR Analyses
# All centers + demographic controls 
sur1 <- prox ~ factor(reltype) + age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious
sur2 <- afrd ~ factor(reltype) + age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious
sur3 <- cent_rept ~ factor(reltype) + age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious
sur4 <- doc_rept ~ factor(reltype) + age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious
sur5 <- doc_know ~ factor(reltype) + age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious
sur6 <- doc_ppl ~ factor(reltype) + age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious
sur7 <- specialty ~ factor(reltype) + age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious
fitsur <- systemfit(list(proxreg = sur1, afrdreg = sur2, centreptreg = sur3, 
                         docreptreg = sur4, docknowreg = sur5, 
                         docpplreg = sur6, specreg = sur7), data=data.tot.iop.leb,
                    method="SUR")
summary(fitsur, residCov = FALSE, equations = FALSE)
round(fitsur$coefficients, 3)

# Block-bootstrapping for SUR
length(unique(data.tot.iop.leb$providercode_unique))
resamp.block <- function(data){
  lookup <- split(1:nrow(data), data$providercode_unique)
  providers <- names(lookup)
  providers.drawn <- sample(providers, length(unique(data$providercode_unique)), replace = TRUE)
  star <- unlist(lookup[providers.drawn])
  samp.blo.data <- data[star,]
}

coefs.sur <- matrix(NA, 70, 1000)
for(i in 1:1000){
  resamp.blo.data <- resamp.block(data.tot.iop.leb)
  fitsur.loop <- systemfit(list(proxreg = sur1, afrdreg = sur2, centreptreg = sur3, 
                                docreptreg = sur4, docknowreg = sur5, 
                                docpplreg = sur6, specreg = sur7), data=resamp.blo.data,
                           method="SUR")
  coefs.sur[,i] <- fitsur.loop$coefficients
  if(i %% 100==0) {
    cat(paste0("iteration: ", i, "\n"))
  }
}

mean.coefs.sur <- apply(coefs.sur, 1, mean, na.rm=T)
sd.coefs.sur <- apply(coefs.sur, 1, sd, na.rm=T)
cbind(mean.coefs.sur, sd.coefs.sur)

# Only sectarian centers, full data
sur1 <- prox ~ outin + age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious
sur2 <- afrd ~ outin + age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious
sur3 <- cent_rept ~ outin + age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious
sur4 <- doc_rept ~ outin + age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious
sur5 <- doc_know ~ outin + age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious
sur6 <- doc_ppl ~ outin + age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious
sur7 <- specialty ~ outin + age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious
fitsur <- systemfit(list(proxreg = sur1, afrdreg = sur2, centreptreg = sur3, 
                         docreptreg = sur4, docknowreg = sur5, 
                         docpplreg = sur6, specreg = sur7), data=data.tot.iops.leb,
                    method="SUR")
summary(fitsur, residCov = FALSE, equations = FALSE)

# Block-bootstrapping for SUR
length(unique(data.tot.iops.leb$providercode_unique))
resamp.block <- function(data){
  lookup <- split(1:nrow(data), data$providercode_unique)
  providers <- names(lookup)
  providers.drawn <- sample(providers, length(unique(data$providercode_unique)), replace = TRUE)
  star <- unlist(lookup[providers.drawn])
  samp.blo.data <- data[star,]
}

coefs.sur <- matrix(NA, 63, 1000)
for(i in 1:1000){
  resamp.blo.data <- resamp.block(data.tot.iops.leb)
  fitsur.loop <- systemfit(list(proxreg = sur1, afrdreg = sur2, centreptreg = sur3, 
                                docreptreg = sur4,  docknowreg = sur5, 
                                docpplreg = sur6, specreg = sur7), data=resamp.blo.data,
                           method="SUR")
  coefs.sur[,i] <- fitsur.loop$coefficients
  if(i %% 100==0) {
    cat(paste0("iteration: ", i, "\n"))
  }
}

mean.coefs.sur <- apply(coefs.sur, 1, mean, na.rm=T)
sd.coefs.sur <- apply(coefs.sur, 1, sd, na.rm=T)
cbind(mean.coefs.sur, sd.coefs.sur)

# With the constrained dataset (patients with access to both ingroup and outgroup centers)
sur1b <- prox ~ outin + age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious
sur2b <- afrd ~ outin + age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious
sur3b <- cent_rept ~ outin + age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious
sur4b <- doc_rept ~ outin + age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious
sur5b <- doc_know ~ outin + age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious
sur6b <- doc_ppl ~ outin + age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious
sur7b <- specialty ~ outin + age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious
fitsur.b <- systemfit(list(proxreg = sur1b, afrdreg = sur2b, centreptreg = sur3b, 
                           docreptreg = sur4b, docknowreg = sur5b, 
                           docpplreg = sur6b, specreg = sur7b), data=data.tot.iops.leb[(data.tot.iops.leb$access_io==1),],
                      method="SUR")

summary(fitsur.b, residCov = FALSE, equations = FALSE)

# Block-bootstrapping for SUR
length(unique(data.tot.iops.leb[(data.tot.iops.leb$access_io==1),]$providercode_unique))
resamp.block.b <- function(data){
  lookup <- split(1:nrow(data), data$providercode_unique)
  providers <- names(lookup)
  providers.drawn <- sample(providers, length(unique(data$providercode_unique)), replace = TRUE)
  star <- unlist(lookup[providers.drawn])
  samp.blo.data <- data[star,]
}

coefs.sur <- matrix(NA, 63, 1000)
for(i in 1:1000){
  resamp.blo.data <- resamp.block.b(data.tot.iops.leb[(data.tot.iops.leb$access_io==1),])
  fitsur.loop <- systemfit(list(proxreg = sur1b, afrdreg = sur2b, centreptreg = sur3b, 
                                docreptreg = sur4b, docknowreg = sur5b, 
                                docpplreg = sur6b, specreg = sur7b), data=resamp.blo.data,
                           method="SUR")
  coefs.sur[,i] <- fitsur.loop$coefficients
  if(i %% 100==0) {
    cat(paste0("iteration: ", i, "\n"))
  }
}

mean.coefs.sur <- apply(coefs.sur, 1, mean, na.rm=T)
sd.coefs.sur <- apply(coefs.sur, 1, sd, na.rm=T)
cbind(mean.coefs.sur, sd.coefs.sur)






#### SECTION: QUALITY DIFFERENCES BETWEEN INGROUP AND OUTGROUP CENTERS ####

## Table 3: Descriptive statistics of the key independent variable and outcome variables in analyses of healthcare quality 
stargazer(data.tot.iops.leb[c("outin", "do_Q12_log", "do_Q14_MC_sum_new","do_Q34_2_log", "do_dei1")], 
          type="text",
          covariate.labels = c("Outgroup dyad", "Questions by the doctor (log)", "Physical examinations", "Min. of examination (log)", 
                               "Doctor effort index"))

## Online Appendix Section A5: Testing the validity of the objective quality variable
# Table A5.1: Is doctor effort index a valid way to measure objective quality?  Association with health quality indicators
fit.e <- lm(data=data.tot.iops.leb, do_dei1 ~ pei_genhealth)
fit.f <- lm(data=data.tot.iops.leb, do_dei1 ~ pei_genhealth + factor(do_Q8_17))
fit.g <- lm(data=data.tot.iops.leb, do_dei1 ~ pei_todayhealth)
fit.h <- lm(data=data.tot.iops.leb, do_dei1 ~ pei_todayhealth + factor(do_Q8_17))
fit.i <- lm(data=data.tot.iops.leb, do_dei1 ~ do_symptom_count_2)
fit.j <- lm(data=data.tot.iops.leb, do_dei1 ~ do_symptom_count_2 + factor(do_Q8_17))
fit.k <- lm(data=data.tot.iops.leb, do_dei1 ~ factor(visit_type_2))
fit.l <- lm(data=data.tot.iops.leb, do_dei1 ~ factor(visit_type_2) + factor(do_Q8_17))

se.e <- list(coef_test(fit.e, vcov = "CR2", 
                       cluster = data.tot.iops.leb$providercode_unique, test = "Satterthwaite")[,2])
se.f <- list(coef_test(fit.f, vcov = "CR2", 
                       cluster = data.tot.iops.leb$providercode_unique, test = "Satterthwaite")[,2])
se.g <- list(coef_test(fit.g, vcov = "CR2", 
                       cluster = data.tot.iops.leb$providercode_unique, test = "Satterthwaite")[,2])
se.h <- list(coef_test(fit.h, vcov = "CR2", 
                       cluster = data.tot.iops.leb$providercode_unique, test = "Satterthwaite")[,2])
se.i <- list(coef_test(fit.i, vcov = "CR2", 
                       cluster = data.tot.iops.leb$providercode_unique, test = "Satterthwaite")[,2])
se.j <- list(coef_test(fit.j, vcov = "CR2", 
                       cluster = data.tot.iops.leb$providercode_unique, test = "Satterthwaite")[,2])
se.k <- list(coef_test(fit.k, vcov = "CR2", 
                       cluster = data.tot.iops.leb$providercode_unique, test = "Satterthwaite")[,2])
se.l <- list(coef_test(fit.l, vcov = "CR2", 
                       cluster = data.tot.iops.leb$providercode_unique, test = "Satterthwaite")[,2])

stargazer(fit.e, fit.f, fit.g, fit.h, fit.i, fit.j, fit.k, fit.l, 
          se = c(se.e, se.f, se.g, se.h, se.i, se.j, se.k, se.l), 
          type="text",
          dep.var.labels = "Doctor effort index",
          covariate.labels = c("General health", "Today's health", "Symptom count", "Visit type: Other", "Visit type: Pregnancy", "Visit type: First-time primary care", "Visit type: Vaccination"),
          omit = "do_Q8_17",
          add.lines =  list(c("Enumerator F.E.", "No", "Yes", "No", "Yes", "No", "Yes", "No", "Yes")),
          notes.align = "l",
          notes.append = FALSE,
          notes = c("* denotes p-value < 0.1, ** denotes p-value < 0.05 and *** denotes p-value < 0.01.", 
                    "Cluster-robust standard errors (at the provider level) are presented in parentheses."))



# Table A5.2: Is doctor effort index a valid way to measure objective quality?  Association with health quality indicators
fit.a <- lm(data=data.tot.iops.leb, pei_subjsat1 ~ do_dei1)
fit.b <- lm(data=data.tot.iops.leb, pei_subjsat1 ~ do_dei1 + age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious +
              factor(visit_type_2))
fit.c <- lm(data=data.tot.iops.leb, pei_subjsat1 ~ do_dei1 + age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious +
              factor(visit_type_2) +
              prox + afrd + doc_rept + cent_rept + doc_know + doc_ppl + specialty)

se.a <- list(coef_test(fit.a, vcov = "CR2", 
                       cluster = data.tot.iops.leb$providercode_unique, test = "Satterthwaite")[,2])
se.b <- list(coef_test(fit.b, vcov = "CR2", 
                       cluster = data.tot.iops.leb$providercode_unique, test = "Satterthwaite")[,2])
se.c <- list(coef_test(fit.c, vcov = "CR2", 
                       cluster = data.tot.iops.leb$providercode_unique, test = "Satterthwaite")[,2])

stargazer(fit.a, fit.b, fit.c, 
          type = "text",
          covariate.labels = c("Doctor effort index", "Patient age", "Respondent age", "Patient gender: Female", "Respondent gender: Female", 
                               "Socioecon. status", "Patient general health", "Religiosity", "Visit type: Other", "Visit type: Pregnancy", "Visit type: Primary", "Visit type: Vaccination",
                               "Proximity", "Affordability", "Doc reputation", "Center reputation", 
                               "Knowing medical staff", "Knowing brokers", "Special care"),
          dep.var.labels = c("Subjective patient satisfaction index"),
          notes.align = "l",
          notes.append = FALSE,
          notes = c("* denotes p-value < 0.1, ** denotes p-value < 0.05 and *** denotes p-value < 0.01.", 
                    "Cluster-robust standard errors (at the provider level) are presented in parentheses."))


jtools::effect_plot(fit.a, pred="do_dei1",
                    x.label = "Doctor effort index",
                    y.label = "Subjective satisfaction index",
                    plot.points = T,
                    interval = T, 
                    cluster = "providercode_unique",
                    robust = "HC2",
                    colors = "indianred") + ylim(-4,1)
jtools::effect_plot(fit.b, pred="do_dei1",
                    x.label = "Doctor effort index",
                    y.label = "Subjective satisfaction index",
                    plot.points = T,
                    interval = T, 
                    cluster = "providercode_unique",
                    robust = "HC2",
                    colors = "blue")+ ylim(-4,1)
jtools::effect_plot(fit.c, pred="do_dei1",
                    x.label = "Doctor effort index",
                    y.label = "Subjective satisfaction index",
                    plot.points = T,
                    interval = T, 
                    cluster = "providercode_unique",
                    robust = "HC2",
                    colors = "green")+ ylim(-4,1)

## Online Appendix Section 6: Control variables
# Control variables: Demographic variables
stargazer(data.tot.iops.leb[c("age_patient", "age_respondent","pei_gender_patient", "pei_gender_respondent", 
                              "ses4", "pei_genhealth", "religious")], 
          type="text")

t.test(data.tot.iops.leb$religious[data.tot.iops.leb$outin==1], 
       data.tot.iops.leb$religious[data.tot.iops.leb$outin==0])

# Control variables: Visit type
table(data.tot.iops.leb$inout, data.tot.iops.leb$visit_type_2, useNA = "always")
round(prop.table(table(data.tot.iops.leb$inout, data.tot.iops.leb$visit_type_2), margin=1), digits=3) 
chisq.test(data.tot.iops.leb$inout, data.tot.iops.leb$visit_type_2, correct=F)

# Control variables: Reasons of selection variables
stargazer(data.tot.iops.leb[c("prox", "afrd","doc_rept", "cent_rept", 
                              "doc_know", "doc_ppl", "specialty", "knoworg_type_loose")], 
          type="text",
          covariate.labels = c("Proximity", "Affordability", "Doc reputation", "Center reputation", 
                               "Knowing medical staff", "Knowing brokers", "Special care", "Knowing the organization"))

## Table 4: Is there an outgroup disadvantage in objective quality? 
fit.obj1 <- lm(do_Q12_log ~ outin +
                 age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious +
                 factor(visit_type_2) + 
                 prox + afrd + doc_rept + cent_rept + doc_know + doc_ppl + specialty + knoworg_type_loose + 
                 factor(do_enum_factor),
               data=data.tot.iops.leb)
fit.obj2 <- lm(do_Q14_MC_sum_new ~ outin + 
                 age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious +
                 factor(visit_type_2) + 
                 prox + afrd + doc_rept + cent_rept + doc_know + doc_ppl + specialty +knoworg_type_loose + 
                 factor(do_enum_factor),
               data=data.tot.iops.leb) 
fit.obj3 <- lm(do_Q34_2_log ~ outin + 
                 age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious +
                 factor(visit_type_2) + 
                 prox + afrd + doc_rept + cent_rept + doc_know + doc_ppl + specialty +knoworg_type_loose + 
                 factor(do_enum_factor),
               data=data.tot.iops.leb)
fit.obj4 <- lm(do_dei1 ~ outin + 
                 age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious +
                 factor(visit_type_2) + 
                 prox + afrd + doc_rept + cent_rept + doc_know + doc_ppl + specialty +knoworg_type_loose + 
                 factor(do_enum_factor),
               data=data.tot.iops.leb)


se.obj1 <- list(coef_test(fit.obj1, vcov = "CR2", 
                          cluster = data.tot.iops.leb$providercode_unique, test = "Satterthwaite")[,2])
se.obj2 <- list(coef_test(fit.obj2, vcov = "CR2", 
                          cluster = data.tot.iops.leb$providercode_unique, test = "Satterthwaite")[,2])
se.obj3 <- list(coef_test(fit.obj3, vcov = "CR2", 
                          cluster = data.tot.iops.leb$providercode_unique, test = "Satterthwaite")[,2])
se.obj4 <- list(coef_test(fit.obj4, vcov = "CR2", 
                          cluster = data.tot.iops.leb$providercode_unique, test = "Satterthwaite")[,2])


stargazer(fit.obj1, fit.obj2, fit.obj3, fit.obj4,
          se = c(se.obj1, se.obj2, se.obj3, se.obj4), 
          type="text",
          covariate.labels = c("Outgroup", 
                               "Patient age", "Respondent age", "Patient gender: Female", "Respondent gender: Female", "Socioecon. status", "Patient general health", "Religiosity", 
                               "Visit type: Other", "Visit type: Pregnancy", "Visit type: Primary", "Visit type: Vaccination",
                               "Proximity", "Affordability", "Doc reputation", "Center reputation", 
                               "Knowing medical staff", "Knowing brokers", "Special care", "Knowing the organization"),
          add.lines =  list(c("Enumerator f.e.", "Yes", "Yes", "Yes", "Yes")),
          omit = "do_enum_factor",
          dep.var.labels = c("Questions asked by doctor (log)", "Physical examinations", "Minutes of exam (log)", "Doctor effort index"),
          notes.align = "l",
          notes.append = FALSE,
          notes = c("* denotes p-value < 0.1, ** denotes p-value < 0.05 and *** denotes p-value < 0.01.", 
                    "Cluster-robust standard errors (at the provider level) are presented in parentheses."))


## Online Appendix A7: Results exploring objective quality differences with weighted data
# Aggregate outgroup variable -- with center selection factor variables -- weighted [Online Appendix]
fit.obj1 <- lm(do_Q12_log ~ outin +
                 age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious +
                 factor(visit_type_2) + 
                 prox + afrd + doc_rept + cent_rept + doc_know + doc_ppl + specialty + knoworg_type_loose + 
                 factor(do_enum_factor),
               data=data.tot.iops.leb, weights = wt)
fit.obj2 <- lm(do_Q14_MC_sum_new ~ outin + 
                 age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious +
                 factor(visit_type_2) + 
                 prox + afrd + doc_rept + cent_rept + doc_know + doc_ppl + specialty + knoworg_type_loose + 
                 factor(do_enum_factor),
               data=data.tot.iops.leb, weights = wt) 
fit.obj3 <- lm(do_Q34_2_log ~ outin + 
                 age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious +
                 factor(visit_type_2) + 
                 prox + afrd + doc_rept + cent_rept + doc_know + doc_ppl + specialty + knoworg_type_loose + 
                 factor(do_enum_factor),
               data=data.tot.iops.leb, weights = wt)
fit.obj4 <- lm(do_dei1 ~ outin + 
                 age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious +
                 factor(visit_type_2) + 
                 prox + afrd + doc_rept + cent_rept + doc_know + doc_ppl + specialty + knoworg_type_loose + 
                 factor(do_enum_factor),
               data=data.tot.iops.leb, weights = wt)


se.obj1 <- list(coef_test(fit.obj1, vcov = "CR2", 
                          cluster = data.tot.iops.leb$providercode_unique, test = "Satterthwaite")[,2])
se.obj2 <- list(coef_test(fit.obj2, vcov = "CR2", 
                          cluster = data.tot.iops.leb$providercode_unique, test = "Satterthwaite")[,2])
se.obj3 <- list(coef_test(fit.obj3, vcov = "CR2", 
                          cluster = data.tot.iops.leb$providercode_unique, test = "Satterthwaite")[,2])
se.obj4 <- list(coef_test(fit.obj4, vcov = "CR2", 
                          cluster = data.tot.iops.leb$providercode_unique, test = "Satterthwaite")[,2])


stargazer(fit.obj1, fit.obj2, fit.obj3, fit.obj4,
          se = c(se.obj1, se.obj2, se.obj3, se.obj4), 
          type="text",
          covariate.labels = c("Outgroup", 
                               "Patient age", "Respondent age", "Patient gender: Female", "Respondent gender: Female", "Socioecon. status", "Patient general health", "Religiosity", 
                               "Visit type: Other", "Visit type: Pregnancy", "Visit type: Primary", "Visit type: Vaccination",
                               "Proximity", "Affordability", "Doc reputation", "Center reputation", 
                               "Knowing medical staff", "Knowing brokers", "Special care", "Knowing the organization"),
          add.lines =  list(c("Enumerator f.e.", "Yes", "Yes", "Yes", "Yes")),
          omit = "do_enum_factor",
          dep.var.labels = c("Questions asked by doctor (log)", "Physical examinations", "Minutes of exam (log)", "Doctor effort index"),
          notes.align = "l",
          notes.append = FALSE,
          notes = c("* denotes p-value < 0.1, ** denotes p-value < 0.05 and *** denotes p-value < 0.01.", 
                    "Cluster-robust standard errors (at the provider level) are presented in parentheses.",
                    "Models are estimated with weights."))


## Online Appendix A8: Using the know-do gap 
data.tot.iops.leb$doc_knowdogap <- as.numeric(data.tot.iops.leb$doc_knowdogap)
stargazer(data.tot.iops.leb[c("doc_knowdogap")], type="text")

fit.obj1 <- lm(doc_knowdogap ~ outin + 
                 age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious +
                 factor(visit_type_2) + 
                 factor(do_enum_factor),
               data=data.tot.iops.leb)
fit.obj2 <- lm(doc_knowdogap ~ outin + 
                 age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious +
                 factor(visit_type_2) + 
                 factor(do_enum_factor),
               data=data.tot.iops.leb, weights = data.tot.iops.leb$wt)

se.obj1 <- list(coef_test(fit.obj1, vcov = "CR2", 
                          cluster = data.tot.iops.leb$providercode_unique, test = "Satterthwaite")[,2])
se.obj2 <- list(coef_test(fit.obj2, vcov = "CR2", 
                          cluster = data.tot.iops.leb$providercode_unique, test = "Satterthwaite")[,2])

stargazer(fit.obj1, fit.obj2, 
          se = c(se.obj1, se.obj2), 
          type="text",
          covariate.labels = c("Outgroup", 
                               "Patient age", "Respondent age", "Patient gender: Female", "Respondent gender: Female", "Socioecon. status", "Patient general health", "Religiosity", 
                               "Visit type: Other", "Visit type: Pregnancy", "Visit type: Primary", "Visit type: Vaccination"),
          add.lines =  list(c("Enumerator f.e.", "Yes", "Yes")),
          omit = "do_enum_factor",
          dep.var.labels = c("Know-do gap"),
          notes.align = "l",
          notes.append = FALSE,
          notes = c("* denotes p-value < 0.1, ** denotes p-value < 0.05 and *** denotes p-value < 0.01.", 
                    "Cluster-robust standard errors (at the provider level) are presented in parentheses.",
                    "Model 2 is estimated with survey design weights that correct for the underrepresentation of party-affiliated centers."))

#### DISCUSSION ####

## Table A9.1
# Variables to test potential mechanisms
stargazer(data.tot.iops.leb[c("do_Q13_log", "socialnetwork2", "pei_docaltruism")], 
          type="text",
          covariate.labels = c("Questions by the patient (log)", "Shared social network mentioned",  "Patient opinion: Altruistic doc"))
t.test(data.tot.iops.leb$pei_docaltruism[data.tot.iops.leb$outin==1], 
       data.tot.iops.leb$pei_docaltruism[data.tot.iops.leb$outin==0])

## For Table A9.2
# Run separately for each outcome variable and for each potential mediator 
fit.deneme <- lm(do_dei1 ~ outin  + factor(do_enum_factor) +
                   age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious +
                   factor(visit_type_2) + 
                   prox + afrd + doc_rept + cent_rept + doc_know + doc_ppl + specialty + knoworg_type_loose,
                 data=data.tot.iops.leb)
med.fit.e <- lm(socialnetwork2 ~  outin  + factor(do_enum_factor) +
                  age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious +
                  factor(visit_type_2) + 
                  prox + afrd + doc_rept + cent_rept + doc_know + doc_ppl + specialty+ knoworg_type_loose,
                data=data.tot.iops.leb[(is.na(data.tot.iops.leb$do_dei1)==0 & is.na(data.tot.iops.leb$socialnetwork2)==0),])
out.fit.e <-  lm(do_dei1 ~ outin  + factor(do_enum_factor) + socialnetwork2 +
                   age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious +
                   factor(visit_type_2) +
                   prox + afrd + doc_rept + cent_rept + doc_know + doc_ppl + specialty + knoworg_type_loose,
                 data=data.tot.iops.leb)

stargazer(fit.deneme, med.fit.e, out.fit.e, type="text", omit = "do_enum_factor")
med.out.e <- mediate(med.fit.e, out.fit.e, treat = "outin", mediator = "socialnetwork2", sims=1000, 
                     conf.level=0.95)
summary(med.out.e)

## For Figure A9.1
# Sensitivity analysis
sens.out.e <- medsens(med.out.e, rho.by=0.1, effect.type="indirect", sims=1000)
summary(sens.out.e)
plot(sens.out.e, main="Sensitivity with Respect to Error Correlation")

## Table A9.3
# Descriptive statistics
stargazer(data.tot.iops.syr[c("outin")], 
          type = "text",
          covariate.labels = c("Outgroup dyad"))

## Table A9.4a
# Objective quality indicators
stargazer(data.tot.iops.syr[c("do_Q12_log", "do_Q14_MC_sum_new","do_Q34_2_log", "do_dei1")], 
          type = "text",
          covariate.labels = c("Questions by the doctor (log)", "Physical examinations", "Min. of examination (log)", "Doctor effort index"))


## Table A9.4b
# Control variables: Demographic variables
stargazer(data.tot.iops.syr[c("age_patient", "age_respondent","pei_gender_patient", "pei_gender_respondent", 
                              "ses4", "pei_genhealth", "religious")], 
          type="text",
          covariate.labels = c("Patient age", "Respondent age", "Patient gender: Female", "Respondent gender: Female", 
                               "Socioecon. status", "Patient general health", "Religiosity"))



## Table A9.4c
# Control variables: Reasons of selection variables
stargazer(data.tot.iops.syr[c("prox", "afrd","doc_rept", "cent_rept", 
                              "doc_know", "doc_ppl", "specialty", "knoworg_type_loose")], 
          type="text",
          covariate.labels = c("Proximity", "Affordability", "Doc reputation", "Center reputation", 
                               "Knowing medical staff", "Knowing brokers", "Special care", "Knowing organization"))

## Table A9.5
# Objective quality with aggregated outgroup variable: Syrian patients
fit.obj1 <- lm(do_Q12_log ~ outin + 
                 age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious +
                 factor(visit_type_2) + 
                 prox + afrd + doc_rept + cent_rept + doc_know + doc_ppl + specialty + knoworg_type_loose + 
                 factor(do_Q8_17),
               data=data.tot.iops.syr)
fit.obj2 <- lm(do_Q14_MC_sum_new ~ outin + 
                 age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious +
                 factor(visit_type_2) + 
                 prox + afrd + doc_rept + cent_rept + doc_know + doc_ppl + specialty + knoworg_type_loose + 
                 factor(do_Q8_17),
               data=data.tot.iops.syr) 
fit.obj3 <- lm(do_Q34_2_log ~ outin + 
                 age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious +
                 factor(visit_type_2) + 
                 prox + afrd + doc_rept + cent_rept + doc_know + doc_ppl + specialty + knoworg_type_loose + 
                 factor(do_Q8_17),
               data=data.tot.iops.syr)
fit.obj4 <- lm(do_dei1 ~ outin + 
                 age_patient + age_respondent + pei_gender_patient + pei_gender_respondent + ses4 + pei_genhealth + religious +
                 factor(visit_type_2) + 
                 prox + afrd + doc_rept + cent_rept + doc_know + doc_ppl + specialty + knoworg_type_loose + 
                 factor(do_Q8_17),
               data=data.tot.iops.syr)


se.obj1 <- list(coef_test(fit.obj1, vcov = "CR2", 
                          cluster = data.tot.iops.syr$providercode_unique, test = "Satterthwaite")[,2])
se.obj2 <- list(coef_test(fit.obj2, vcov = "CR2", 
                          cluster = data.tot.iops.syr$providercode_unique, test = "Satterthwaite")[,2])
se.obj3 <- list(coef_test(fit.obj3, vcov = "CR2", 
                          cluster = data.tot.iops.syr$providercode_unique, test = "Satterthwaite")[,2])
se.obj4 <- list(coef_test(fit.obj4, vcov = "CR2", 
                          cluster = data.tot.iops.syr$providercode_unique, test = "Satterthwaite")[,2])


stargazer(fit.obj1, fit.obj2, fit.obj3, fit.obj4,
          se = c(se.obj1, se.obj2, se.obj3, se.obj4), 
          type="text",
          covariate.labels = c("Outgroup", 
                               "Patient age", "Respondent age", "Patient gender: Female", "Respondent gender: Female", "Socioecon. status", "Patient general health", "Religiosity", 
                               "Visit type: Other", "Visit type: Pregnancy", "Visit type: Primary", "Visit type: Vaccination",
                               "Proximity", "Affordability", "Doc reputation", "Center reputation", 
                               "Knowing medical staff", "Knowing brokers", "Special care", "Knowing the organization"),
          add.lines =  list(c("Enumerator f.e.", "Yes", "Yes", "Yes", "Yes")),
          omit = "do_Q8_17",
          dep.var.labels = c("Questions asked by doctor (log)", "Physical examinations", "Minutes of exam (log)", "Doctor effort index"),
          notes.align = "l",
          notes.append = FALSE,
          notes = c("* denotes p-value < 0.1, ** denotes p-value < 0.05 and *** denotes p-value < 0.01.", 
                    "Cluster-robust standard errors (at the provider level) are presented in parentheses."))









